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1  Introduction 


The  extraction  of  molecular  Hamiltonian  information  from  laboratory  data  has  been 
a  goal  of  considerable  importance  in  chemistry.  In  general,  the  data  can  come  from 
either  stationary  or  time-dependent  experiments,  with  the  latter  just  recently  becom¬ 
ing  available  at  interesting  ultrashort  time  scales[l].  In  addition,  there  are  numerous 
spectroscopic,  and  other  sources,  of  time-independent  data.  A  significant  bottleneck 
in  this  field  has  been  the  lack  of  reliable  algorithms  for  inverting  such  data  to  their 
underlying  physical  (Hamiltonian)  information[2,  3,  4,  5,  6]. 

The  algorithm  presented  in  this  work  may  be  traced  back  to  tracking  control  of 
molecular  motion[7].  In  the  inversion  context,  a  ’’track”  <  0{t)  >  was  assumed 
to  be  observed  laboratory  data,  and  the  goal  was  to  find  the  potential  or  electric 
dipole  function  consistent  with  the  track.  This  report  presents  the  basic  rudiments 
of  the  idea  that  was  developed  under  this  project,  and  the  results  suggested  that  the 
concepts  were  quite  well-founded. 

V 

A  variety  of  dynamical  experiments  [8]  and  theoretical  techniques  provided  the 
basis  for  the  inversion  research.  The  last  three  decades  witnes^d  the  development  of 
time-dependent  quantum-mechanical  methods  for  molecular  dynamics  [9],  which  were 
made  possible  by  the  efficient  numerical  algorithms  of  wavepacket  propagation[10]. 
Combining  the  advances  in  both  the  experimental  and  the  numerical  techniques,  we 
proposed  a  new  inversion  algorithm  that  reconstructs  the  underlying  potential  energy 
surfaces  (PES)  and  dipole  function  from  temporal  measurements. 


2  Methodology 

Different  types  of  observables  and  different  initial  conditions  can  arise,  depending  on 
the  nature  of  the  experiments.  Here,  we  consider  a  single  observable  operator  0  with 
a  sequence  of  initial  conditions  ipj{t  =  0),j  =  1,2. ..No-  The  Hermitian  operator 
0  corresponded  to  a  physical  quantity  observed  in  laboratory  as  a  function  of  time. 
For  example,  in  the  case  of  position  tracking,  0  =  x,  and  <  x{t)  >  is  the  measured 
mean  position  of  the  atoms  in  the  molecule  (e.g.,  measurement  of  ultrafast  X-ray  or 
electron  diffraction). 

The  experiments  were  assumed  to  provide  measurements  <  Oj{t)  >=<  Tpj{t)\0\ipj 
(t)  >,  j  =  1,  2...Nd-  According  to  the  Heisenberg  equation  of  motion,  we  had 

=  <  V]\ijj  >  -€j{t)  <  >  (1) 

where  H  =  T  +  V  —  was  the  Hamiltonian  for  the  system,  T  was  the  kinetic 
energy  operator,  /x(x)  was  the  dipole  moment  and  e{t)  was  an  external  field. 

For  any  arbitrary  Hermitian  operator  0,  Eq.  (1)  could  be  iecast  into  the  form  of 
a  Fredholm  integral  equation  of  the  first  kind  for  the  potential 

j  K{j,x)V[x)dx  =  gj  (2) 


where 


9j  =  >, 


K{j,x)  = 


A  similar  expression  could  be  written  for  the  dipole  function,  if  it  was  also  sought. 
The  kernels  K{j,x)  in  Eq.  (3)  implicitly  depended  on  time.  We  assumed  that  [O,  V\ 
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was  not  identically  zero.  If  [0,  V]  =  0,  then  Eq.  (1)  was  not  an  invertible  relation 
as  it  stood.  In  such  a  case,  the  differentiation  of  the  track  <  Oj{t)  >  with  respect 
to  time  must  be  continued  to  obtain  a  higher-order  equation  which  is  invertible[ll]. 
An  example  was  the  case  of  position  tracking  0  —  x.  Although  [x,  V(x)]  =  0,  the 
invertible  relation  for  V  was  obtained  after  differentiating  <  x{t)  >  twice  with  respect 
to  time. 

Given  the  above  basic  relations,  a  family  of  algorithms  was  suggested  to  achieve 
the  desired  inversion.  Many  issues  arose  in  this  regard,  and  only  the  simplest  approach 
was  explored  in  the  research.  The  procedure  for  the  potential  inversion  is  given  below 
in  detail,  based  on  the  field  s(t)  being  zero.  This  discussion  can  be  easily  adapted 
to  the  dipole,  inversion  case.  The  actual  laboratory  data  <  Oj{t)  >  is  a  convolution  -- 
of  the  wavefunction  in  configuration  space  (or  some  Other  representation),  and  thus, 
no  urdque  wavefunction  ipj  corresponded  to  the  data.  To  start  the  inversion,  a  set  of 
initial  conditions  =  0)  was  needed  to  be  consistent  with  <  ‘ipj{0)\O\‘tjjj{0)  >=< 
Oj{t  =  0)  >,  and  much  freedom  existed  in  the  choice.  With  the  known  initial  states 
=  0)  and  the  initial  data  <  Oj{t  =  0)  >,  the  unknown  potential  V  was  the 
solution  of  the  integral  equation  (2).  Since  the  initial  nonstationary  states  ^j(O)  are 
usually  localized  in  the  Franck- Condon  region  in  an  electronic  excitation  process,  the 
solution  of  the  integral  equation  consequently  only  gave  the  local  potential  Vq  in  the  ' 
Franck-Condon  region.  However,  this  piece  of  potential  Vq  could  be  used  to  propagate 
the  wavefunctions  ipj{0)  to  the  next  time  fpjiAt)  according  to  the  time-dependent 
Schrodinger  equation, 

=  (T +  V)\i,i  >  -■  (4) 

Then,  the  solution  of  the  integral  equation  (2)  gave  the  next  piece  of  the  potential 
With  V^t,  the  wavefunctions  at  the  next  time  step  V’j(2Ai)  could  obtained  by 
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propagation,  which  allows  for  solving  the  integral  equation  to  get  V2Ati  etc.  This 
procedure  could  be  effected  in  terms  of  the  split  operator  method[12].  The  process  of 
inversion  plus  propagation  continued  until  the  wavefunctions  sampled  all  dynamically 
accessible  regions  of  the  PES.  Finally,  the  potential  V  was  the  union  of  the  sampled 
pieces  along  the  track. 

3  The  Need  for  Regularization 

The  solution  of  the  integral  equation  (2)  was  generally  not  unique.  This  was  due  to 
the  limited  information  we  could  get  from  the  inevitably  incomplete  data.  In  addition, 
the  solution  was  likely  not  stable  against  small  changes  in  the  data,  which  were  known 
as  ill-posedness.  Both  of  these  issues  are  common  to  virtually  all  inverse  problems. 
However,  we  did  not  seek  an  arbitrary  potential  or  dipole  function,  but  a  physical  one 
which  was  capable  of  reproducing  the  data  and  also  met  some  fundamental  physical 
requirements  and  boundary  conditiqjis.  These  requirements  and  conditions  consti¬ 
tuted  our  prior  knowledge  about  the  potential  or  dipole  function.  For  example,  the 
physical  PES  should  be  a  smooth  function  and  negligible  asjmptoticaUy  as  a;  ^  oo. 
As  a  general  inverse-problem  approach  was  adopted  in  this  research,  we  did  not  as¬ 
sume  any  biased  functional  form  for  the  potential.  Without  further  restriction,  this  ^ 
allowed  multiple  solutions  satisfying  Eq.  (2),  with  some  of  them  highly  irregular  and 
even  singular.  The  incorporation  of  prior  knowledge  provided  some  mild  general  con¬ 
straints  on  the  solutions.  Thus,  by  incorporating  this  prior  knowledge,  the  solution 

could  be  stabihzed,  which  is  referred  as  a  regularization  procedure.  Here,  we  used  the 

i 

very  mild  criterion  of  smoothness  of  the  potential  as  our  prior  knowledge  to  regularize 
the  solution. 
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In  order  to  implement  the  regularization,  Eq.  (2)  was  rewritten  as  [2] 


J  {j,  {x)dx  =  gj 


(5) 


where 


= 


iPV 
dx^  ’ 


K^^'^{j,x)  =  -  [  ^\j,x')dx'  (6) 

The  regularized  solution  of  Eq.  (2)  was  achieved  by  minimizing  the  following  functional[2] 


J 


leading  to 


=  II  / isrWy, x)V»(x)dx  -  gM  +  a||VW(i)|| 


n\ 


(7) 


(8) 


where  Uj,Vj{x)  were  the  singular  functions  for  the  kernel  in  Eq.  (5),  with  singular 


values  u)j. 


The  regularization  parameter  a  denoted  the  trade-off  between  reproducing  the 
data  and  the  solution  being  smooth.  We  constructed  an  error  function  of  a, 


Error{a)  =  \\  j x)V^'"^{x)dx  -  gj\\ 


(9) 


The  optimal  a*  was  located  at  the  global  minimum  of  this  error  function. 

The  inversion  was  actually  achieved  through  the  derivative  which 

imposed  a  stronger  smoothness  condition  on  V{x)  for  n  >  1.  To  obtain  V{x),  we 
needed  to  integrate  n  times 


1*00 

^  _  /  dx'V^^^ix')  - 
Jx 


(10) 


The  above  proposed  algorithm  was  presented  in  one  dimension,  but  it  could  be  gen¬ 
eralized  to  higher  dimensions  as  well. 
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4  Example 

We  simulated  the  inversion  of  the  potential  and  dipole  function  of  a  diatomic  molecule. 
Suppose  that  the  mean  positions  <  Xj{t)  >=<  •0(t)|x|^(t)  >  of  the  wavepackets  could 
be  measured  in  the  experiment  (e.g.,  by  ultrafast  X-ray  or  electron  diffraction).  Since 
the  Hesienberg  Eq.  (1)  for  6  =  re  was  not  invertible,  the  second  derivative  could  be 
invoked  to  give 

^  <  Xj{t)  >=  -y  dx\ipj\'^V^'^^{x)  +  j da:e(t)|'0jf /’(a;)  (H) 

where  the  superscript  1  denoted  the  first  order  space  derivative. 

We  proceeded  to  reconstruct  the  potential  and  the  permanent  dipole  functions 
from  the  time-dependent  data  <  Xj  (t)  >  by  solving  the  Fredholm  integral  equation  of 
the  first  kind,  Eq.  (2).  We  applied  the  direct  inversion  procedure  to  a  model  problem 
of  a  bound  diatomic  molecule. 

In  the  simulation,  the  initial  wavefunctions  were  taken  as  Gaussian  wavepackets, 

» 

^qoj,Poj  i^)  =  C'errp  (^^PojX  -  ,  j  =  1, 2,  ...No  (12) 

where  C  =  was  the  normalization  constant.  The  initial  condition  in  Eq. 

(12)  was  consistent  with  the  track  <  Xj{t)  >. 

The  results  presented  below  were  based  on  the  following  conditions.  The  second 
derivative  of  the  potential  was  required  to  be  smooth,  i.e.,  n  =  2.  The  number  of 
data  is  Nd  =  16,  poj  =  0,  qoj  =  —0.9  +  {j  —  1)  Aa;  where  Ax  =  0.01  was  the  numerical 
grid  spacing;  and  Xmin  —  -2,  Xmax  =  8  specified  the  size  of  the  numerical  grid. 

The  inverted  potential  was  compared  with  the  exact  one  V  in  Fig.  la.  The  over¬ 
all  quality  of  the  inverted  potential  was  excellent.  As  a  more  sensitive  verification, 
the  wavefunction  at  fixed  interval  of  time  was  compared  to  those  obtained  by  the 
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propagation  on  the  exax^t  potential  V{x)  (so-called  exact  wavefunctions)  in  Fig.  lb. 
The  wavefunctions  propagated  on  the  inverted  potential  had  no  essential  difference 
from  the  exact  ones,  except  that  a  small  fluctuation  was  visible  at  later  times,  which 
resulted  from  the  finite  size  of  the  chosen  numerical  grid.  A  similar  inversion  was 
carried  out  to  extract  the  dipole  function  in  this  model,  and  the  results  for  the  dipole 
function  are  given  in  Fig.  2.  Similar  quahty  excellent  results  were  obtained  for  the 
inversion. 

5  Conclusion 

This  report  summarizes  a  new  direct  inversion  algorithm  under  development  for  recov¬ 
ering  molecular  potentials  and  dipole  functions  from  sequential  time  measurements. 
The  direct  algorithm  exploits  efficient  numerical  methods  of  wavepacket  propagation. 
Time  is  expHcitly  utilized  as  an  important  parameter  in  organizing  the  inversion.  We 
have  demonstrated  this  algorithm  iii  a  simulation  for  recovering  the  potential  and 
dipole  moment  for  a  model  diatomic  molecule. 

The  proposed  inversion  algorithm  is  direct  and  computafionally  efficient.  The 
traditional  time-independent  inversion  schemes  usually  require  many  iterations  before 
convergence.  This  can  be  very  time  consuming  computationally.  In  the  tracking 
method,  the  inverse  problem  is  solved  just  once  at  each  time  step.  Furthermore, 
the  size  of  the  inverse  problem  solved  in  the  direct  approach  can  be  much  smaller 
{Nd  «  10)  than  that  in  the  time-independent  methods  {No  ^  10^  in  one  dimension), 
because  the  potential  sampled  by  the  localized  wavepackets  can  be  inverted  locally 

i 

which  is  in  contrast  to  the  global  potential  sampled  by  all  the  eigenstates  or  the 
continuum  states.  In  the  direct  method  most  of  the  computational  effort  deals  with 
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the  wavepacket  propagation,  which  is  much  easier. 

The  time-dependent  tracking  approach  to  inversion  also  opens  up  the  possibihty 
to  recover  the  potential  from  time-independent  data,  since  the  latter  data  can  be 
synthesized  from  the  time-independent  information.  This  topic,  as  well  as  a  host 
of  other  issues,  needs  to  be  explored  to  fully  exploit  the  capabihties  of  the  direct 
algorithm.  The  results  of  this  research  opens  up  this  field  for  further  exploration  of 
these  new  concepts  of  laboratory  data  inversion. 
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Fig.  1  :  The  results  for  the  potential  inversion.  The  initial  conditions  are;  ?oi  =  -0-9 
and  po  =  0. 

(a)  V{x).  The  squares  are  the  inverted  results,  and  the  solid  line  is  the  exact 

potential.  - 

(b)  Snapshots  of  at  the  successive  intervals  of  time,  where  the  solid  line 

is  recorded  during  the  inversion  and  dotted  line  is  the  exact  result  obtained 
by  the  propagation  on  Vm{x). 
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